####################################################################################################### 
# Replication Script  
# Anti-Muslim Bias in Foreign Policy Attitudes: Experimental Evidence from 13 European Countries
#######################################################################################################

#######################################################################################################
#### packages ####
set.seed(123)
library(rio)
library(MASS)
library(mokken)
library(lme4)
library(nlme)
library(lmerTest)
library(effectsize)
library(marginaleffects)
library(modelsummary)
library(clubSandwich)
library(estimatr)
library(ggplot2)
library(ggpubr)
library(dplyr)
library(psych)
library(dplyr)
library(texreg)
library(tinytable)

#######################################################################################################

#######################################################################################################
#### sessionInfo() ####

# R version 4.4.0 (2024-04-24)
# Platform: aarch64-apple-darwin20
# Running under: macOS Sonoma 14.5

# Matrix products: default
# BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
# LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

#locale:
#  [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

# time zone: America/Chicago
# tzcode source: internal

# attached base packages:
#  [1] stats     graphics  grDevices utils     datasets  methods   base   

# other attached packages:
#  [1] effectsize_0.8.9       tinytable_0.3.0        texreg_1.39.4          psych_2.4.6.26         dplyr_1.1.4           
#  [6] ggpubr_0.6.0           ggplot2_3.5.1          estimatr_1.0.4         clubSandwich_0.5.11    modelsummary_2.1.1    
# [11] marginaleffects_0.21.0 lmerTest_3.1-3         nlme_3.1-165           lme4_1.1-35.5          Matrix_1.7-0          
# [16] mokken_3.1.2           poLCA_1.6.0.1          scatterplot3d_0.3-44   MASS_7.3-61            rio_1.2.0   

# loaded via a namespace (and not attached):
# [1] mnormt_2.1.1        gridExtra_2.3       sandwich_3.1-0      readxl_1.4.3        rlang_1.1.4         magrittr_2.0.3     
# [7] compiler_4.4.0      mgcv_1.9-1          vctrs_0.6.5         pkgconfig_2.0.3     fastmap_1.2.0       backports_1.5.0    
# [13] labeling_0.4.3      utf8_1.2.4          nloptr_2.1.1        purrr_1.0.2         xfun_0.46           broom_1.0.6        
# [19] parallel_4.4.0      R6_2.5.1            tables_0.9.28       parallelly_1.37.1   car_3.1-2           boot_1.3-30        
# [25] lmtest_0.9-40       cellranger_1.1.0    numDeriv_2016.8-1.1 estimability_1.5.1  Rcpp_1.0.13         knitr_1.48         
# [31] future.apply_1.11.2 zoo_1.8-12          parameters_0.22.1   R.utils_2.12.3      splines_4.4.0       tidyselect_1.2.1   
# [37] rstudioapi_0.16.0   abind_1.4-5         codetools_0.2-20    listenv_0.9.1       lattice_0.22-6      tibble_3.2.1       
# [43] withr_3.0.0         bayestestR_0.14.0   coda_0.19-4.1       future_1.33.2       pillar_1.9.0        carData_3.0-5      
# [49] checkmate_2.3.1     insight_0.20.2      generics_0.1.3      munsell_0.5.1       scales_1.3.0        minqa_1.2.7        
# [55] globals_0.16.3      xtable_1.8-4        glue_1.7.0          emmeans_1.10.3      tools_4.4.0         data.table_1.16.0  
# [61] ggsignif_0.6.4      mvtnorm_1.2-5       cowplot_1.1.3       grid_4.4.0          tidyr_1.3.1         datawizard_0.12.2  
# [67] colorspace_2.1-0    performance_0.12.2  Formula_1.2-5       cli_3.6.3           fansi_1.0.6         gtable_0.3.5       
# [73] R.methodsS3_1.8.2   rstatix_0.7.2       digest_0.6.36       farver_2.1.2        htmltools_0.5.8.1   R.oo_1.26.0        
# [79] lifecycle_1.0.4     httr_1.4.7         
#######################################################################################################

#######################################################################################################
#### Data Import and Check ####

data <- rio::import("Muslim_Bias_Data.csv")

## Create outcome index ##
data$index <- (data$persecution_stop+data$ch_govt_stop+data$govt_pressure +data$un_pressure + data$petition)/5

## check descriptives ## 
data %>% group_by(Country) %>% summarize(
  mean_persecution = mean(persecution_stop),
  mean_ch_govt_stop = mean(ch_govt_stop),
  mean_govt_pressure = mean(govt_pressure),
  mean_un_pressure = mean(un_pressure),
  mean_petition = mean(petition),
  mean_index = mean(index)
)

## Re-order treatment factors ## 
data$sb_religion <- factor(data$sb_religion, levels = c("Muslims", "Buddhists", "Christians", "Taoists"))

#######################################################################################################

#######################################################################################################
#### Mokken Scale Analysis and Cronbach of Outcome Index ####

index <- dplyr::select(data,25:29)

store <- coefH(index)
store$H 
store$Hij

aisp(index, search="ga")

monotonicity <- check.monotonicity(index)
summary(monotonicity)
monotonicity_table <- as.data.frame(summary(monotonicity))

## Table A1 ##
tt(monotonicity_table, digits=2) |> print("latex")

## Table A2 @@
pmatrix <- check.pmatrix(index)
summary(pmatrix)
pmatrix_table <- as.data.frame(summary(pmatrix))
tt(pmatrix_table, digits=2) |> print("latex")

## Cronbach's Alpha ## 
items <- dplyr::select(data, persecution_stop, ch_govt_stop, govt_pressure, un_pressure, petition)

psych::alpha(items, na.rm=T, cumulative = FALSE) 

#######################################################################################################

#######################################################################################################
#### Figure 1: Treatment Effects in the Pooled Sample ####

## Random Intercepts Model ##
model_re <- lmer(index ~ as.factor(sb_religion) + (1 | Country), data = data)
summary(model_re)

## Fixed Effects Model
model_fe <- lm(index ~ as.factor(sb_religion) + as.factor(Country), data = data)
summary(model_fe)

## Create Attention Check Success Variable ##
data$attention_pass <- ifelse(data$sb_religion==data$attention_check,1,0)
table(data$attention_pass)

## Subset on Attention Check Success ##
data_attention <- data[data$attention_pass==1, ]

## Re-estimate both RI and FE model ##
model_re_attention <- lmer(index ~ as.factor(sb_religion) + (1 | Country), data = data_attention)
summary(model_re_attention)

model_fe_attention <- lm(index ~ as.factor(sb_religion) + as.factor(Country), data = data_attention)
summary(model_fe_attention)

models <- list(model_re,model_fe, model_re_attention, model_fe_attention)

## Create Figure ## 
plot_index <- modelplot(models, size=0.3, coef_omit = c(5:18), coef_rename=c("Muslims \n (reference)","Buddhists","Christians","Taoists")) + xlim(0,0.55) + geom_vline(xintercept = 0, linetype = 'dashed') + theme_bw() + xlab("ATE/CATE \n Opposition to Persecution (1-7)") +  theme(axis.text.x = element_text(size = 8), axis.text.y = element_text(size = 8), axis.title.x = element_text(size = 9)) + labs(color="")  + scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9", "#009E73") , guide = guide_legend(reverse = FALSE), labels=c("RI","FE","RI:Attention","FE:Attention")) + theme(legend.position="bottom")

plot_index

#ggsave("Figure1_ATEs.png",plot_index, width=6 ,height = 4.5, dpi=600)

cohens_d_results_1 <- effectsize::cohens_d(index~sb_religion, data=data[data$sb_religion %in% c("Muslims","Buddhists"), ])
cohens_d_results_2 <- effectsize::cohens_d(index~sb_religion, data=data[data$sb_religion %in% c("Muslims","Christians"), ])
cohens_d_results_3 <- effectsize::cohens_d(index~sb_religion, data=data[data$sb_religion %in% c("Muslims","Taoists"), ])

hedges_g_results_1 <- effectsize::hedges_g(index~sb_religion, data=data[data$sb_religion %in% c("Muslims","Buddhists"), ])
hedges_g_results_2 <- effectsize::hedges_g(index~sb_religion, data=data[data$sb_religion %in% c("Muslims","Christians"), ])
hedges_g_results_3 <- effectsize::hedges_g(index~sb_religion, data=data[data$sb_religion %in% c("Muslims","Taoists"), ])

## Table A3: Index Models (Figure 1, Main Text) ## 
table1_models <- list(model_re,model_fe, model_re_attention, model_fe_attention)
modelsummary(table1_models, output = "latex", coef_map = , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))

#######################################################################################################

#######################################################################################################
#### Figure 2: Predicted Opposition Across Countries ####

## Run and clean results for linear regression in each country ## 
country_regs <- list()
countries <- unique(data$Country)

for (country in countries) {
  country_data <- subset(data, Country == country)
  model <- lm(index ~ as.factor(sb_religion), data = country_data)
  country_regs[[country]] <- model
}
country_names <- names(country_regs)
sorted_country_names <- sort(country_names)
country_regs <- country_regs[sorted_country_names]

summary(country_regs$CZ)
summary(country_regs$ES)
summary(country_regs$FR)
summary(country_regs$GE)
summary(country_regs$HU)
summary(country_regs$IT)
summary(country_regs$LV)
summary(country_regs$PL)
summary(country_regs$RU)
summary(country_regs$SK)
summary(country_regs$SR)
summary(country_regs$SW)
summary(country_regs$UK)

## Code for Tables A4 and A5 in SI ## 
modelsummary(country_regs[1:7], output = "latex", coef_map = , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))
modelsummary(country_regs[8:13], output = "latex", coef_map = , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))

## Code for Figure 2 Plot ##
CZ <- avg_predictions(country_regs$CZ, by="sb_religion", conf_level = .95)
CZ$country <- "CZ"

SK <- avg_predictions(country_regs$SK, by="sb_religion", conf_level = .95)
SK$country <- "SK"

HU <- avg_predictions(country_regs$HU, by="sb_religion", conf_level = .95)
HU$country <- "HU"

UK <- avg_predictions(country_regs$UK, by="sb_religion", conf_level = .95)
UK$country <- "UK"

PL <- avg_predictions(country_regs$PL, by="sb_religion", conf_level = .95)
PL$country <- "PL"

IT <- avg_predictions(country_regs$IT, by="sb_religion", conf_level = .95)
IT$country <- "IT"

RU <- avg_predictions(country_regs$RU, by="sb_religion", conf_level = .95)
RU$country <- "RU"

ES <- avg_predictions(country_regs$ES, by="sb_religion", conf_level = .95)
ES$country <- "ES"

SW <- avg_predictions(country_regs$SW, by="sb_religion", conf_level = .95)
SW$country <- "SW"

FR <- avg_predictions(country_regs$FR, by="sb_religion", conf_level = .95)
FR$country <- "FR"

GE <- avg_predictions(country_regs$GE, by="sb_religion", conf_level = .95)
GE$country <- "GE"

SR <- avg_predictions(country_regs$SR, by="sb_religion", conf_level = .95)
SR$country <- "SR"

LV <- avg_predictions(country_regs$LV, by="sb_religion", conf_level = .95)
LV$country <- "LV"

combined_data <- rbind(CZ, SK, HU, UK, PL, IT, RU, ES, SW, FR, GE, SR, LV)

country_plot <- ggplot(combined_data, aes(x = country, y = estimate, color = sb_religion)) + ylab("Predicted Opposition to Persecution (1-7)") + xlab("") +theme_bw() + labs(title="") + theme(plot.title = element_text(size = 8, hjust = 0.5)) + ylim(3.5,5.7) +geom_point(size=1.5, position=position_dodge(width=0.5)) + theme(legend.position="bottom") +  theme(legend.title=element_blank()) + scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9", "#009E73")) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high, xmin = country, xmax = country), width=0, position=position_dodge(width=0.5))

country_plot

#ggsave("Figure2_ATEs_Countries.png",country_plot, width=6.5 ,height = 4.5, dpi=600)


## Table A6: Numerical Estimates for Figure 2  ## 
subset_table <- combined_data[, c(11, 2, 3, 8, 9), drop=FALSE] %>% arrange(country)
latex_code <- knitr::kable(subset_table, format = "latex", col.names = c("Country", "Group", "Estimate", "Conf. Low", "Conf. High"), booktabs=TRUE, linesep="", digits = 2) 
cat(latex_code)

## Table A7-A8 Country-Level Results Attention Check Failures Excluded ##

country_regs <- list()
countries <- unique(data_attention$Country)

for (country in countries) {
  country_data <- subset(data_attention, Country == country)
  model <- lm(index ~ as.factor(sb_religion), data = country_data)
  country_regs[[country]] <- model
}
country_names <- names(country_regs)
sorted_country_names <- sort(country_names)
country_regs <- country_regs[sorted_country_names]

summary(country_regs$CZ)
summary(country_regs$ES)
summary(country_regs$FR)
summary(country_regs$GE)
summary(country_regs$HU)
summary(country_regs$IT)
summary(country_regs$LV)
summary(country_regs$PL)
summary(country_regs$RU)
summary(country_regs$SK)
summary(country_regs$SR)
summary(country_regs$SW)
summary(country_regs$UK)

modelsummary(country_regs[1:7], output = "latex", coef_map = , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))
modelsummary(country_regs[8:13], output = "latex", coef_map = , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))


#######################################################################################################

#######################################################################################################
#### Figure 3: Conditional Average Treatment Effects Across Participants’ Attitudes Towards Muslims ####

table(data$Muslim_thermometer)

model_re_moderation <- lmer(index ~ as.factor(sb_religion)*Muslim_thermometer + (1 | Country), data = data)
summary(model_re_moderation)

model_fe_moderation <- lm(index ~ sb_religion*Muslim_thermometer + as.factor(Country), data = data)
summary(model_fe_moderation)

## Robustness check for additional manipulation concerning Muslim thermometer question ##
model_re_moderation_check <- lmer(index ~ as.factor(sb_religion)*Muslim_thermometer*sb_attitudes + (1 | Country), data = data)
summary(model_re_moderation_check)

model_fe_moderation_check <- lm(index ~ as.factor(sb_religion)*Muslim_thermometer*sb_attitudes + as.factor(Country), data = data)
summary(model_fe_moderation_check)

muslim_moderation_models <- list(model_re_moderation,model_fe_moderation, model_re_moderation_check,model_fe_moderation_check)

## Table A9: Muslim Thermometer Moderation Models ##
modelsummary(muslim_moderation_models, output = "latex", coef_map = , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))

## Figure 3: Predicted Opposition Across Attitudes Towards Muslims ##
prediction_data <- plot_predictions(model_fe_moderation, by=c("Muslim_thermometer","sb_religion"), draw = FALSE)

muslim_therm_plot <- ggplot(prediction_data, aes(x = Muslim_thermometer, y = estimate, color = sb_religion, fill = sb_religion)) +
  geom_smooth(method = "lm", se = TRUE, alpha = 0.3) +
  ylab("Predicted Opposition to Persecution (1-7)") + 
  xlab("Attitudes Towards Muslims") +   
  scale_y_continuous(limits = c(3.85, 5.85), breaks = seq(3.8, 5.9, by = 0.4)) +
  theme_bw() + 
  labs(title="", color="") + 
  theme(plot.title = element_text(size = 8, hjust = 0.5), 
        legend.position="bottom", 
        legend.title = element_blank()) +
  scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9", "#009E73")) +
  scale_fill_manual(values = c("gray","gray","gray","gray")) +
  guides(fill = "none")

muslim_therm_plot

#ggsave("Figure3_Predicted_Muslim_Therm.png",muslim_therm_plot, width=6.5 ,height = 5, dpi=600)

## Table A10: Multilevel Models (Pooled Muslims Attitudes at Country Level) ##

data <- data %>% 
  group_by(Country) %>% 
  mutate(Muslim_thermometer_mean = mean(Muslim_thermometer),
         Muslim_thermometer_centered = Muslim_thermometer - Muslim_thermometer_mean) %>%
  ungroup()


table(data$Muslim_thermometer_mean)
model_ml_index<- lme(index ~ Muslim_thermometer_mean*sb_religion + Muslim_thermometer_centered*sb_religion,
                     random = ~ sb_religion | Country, 
                     data = data,
                     method = "REML",
                     control = lmeControl(opt = "optim"))

model_ml_worry <- lme(persecution_stop ~ Muslim_thermometer_mean*sb_religion + Muslim_thermometer_centered*sb_religion, 
                      random = ~ sb_religion | Country, 
                     data = data,
                     method = "REML",
                     control = lmeControl(opt = "optim"))

model_ml_stop <- lme(ch_govt_stop ~ Muslim_thermometer_mean*sb_religion + Muslim_thermometer_centered*sb_religion, 
                     random = ~ sb_religion | Country, 
                     data = data,
                     method = "REML",
                     control = lmeControl(opt = "optim"))

model_ml_country <- lme(govt_pressure ~ Muslim_thermometer_mean*sb_religion + Muslim_thermometer_centered*sb_religion, 
                     random = ~ sb_religion | Country, 
                     data = data,
                     method = "REML",
                     control = lmeControl(opt = "optim"))

model_ml_UN <- lme(un_pressure ~ Muslim_thermometer_mean*sb_religion + Muslim_thermometer_centered*sb_religion, 
                        random = ~ sb_religion | Country, 
                        data = data,
                        method = "REML",
                        control = lmeControl(opt = "optim"))

model_ml_behavior <- lme(petition ~ Muslim_thermometer_mean*sb_religion + Muslim_thermometer_centered*sb_religion, 
                   random = ~ sb_religion | Country, 
                   data = data,
                   method = "REML",
                   control = lmeControl(opt = "optim"))

VarCorr(model_ml_index)
VarCorr(model_ml_worry)
VarCorr(model_ml_stop)
VarCorr(model_ml_country)
VarCorr(model_ml_UN)
VarCorr(model_ml_behavior)

models <- list(model_ml_index,model_ml_worry, model_ml_stop, model_ml_country, model_ml_UN, model_ml_behavior)



modelsummary(models, output = "latex", coef_map = c(
) , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))
#######################################################################################################

#######################################################################################################
#### Figure 4: Conditional Average Treatment Effects Across Participants’ Attitudes Towards China ####

table(data$China_therm)

model_re_moderation <- lmer(index ~ as.factor(sb_religion)*China_therm + (1 | Country), data = data)
summary(model_re_moderation)

model_fe_moderation <- lm(index ~ sb_religion*China_therm + as.factor(Country), data = data)
summary(model_fe_moderation)

## Table A11: China Attitudes Moderation Models ## 
moderation_models <- list(model_re_moderation, model_fe_moderation)

modelsummary(moderation_models, output = "latex", coef_map = , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('$^a$' = .1, '*'=0.05))


## Figure 4: Predicted Opposition Across Attitudes Towards China ##

prediction_data <- plot_predictions(model_fe_moderation, by=c("China_therm","sb_religion"), draw=FALSE)

china_moderation_plot <- ggplot(prediction_data, aes(x = China_therm, y = estimate, color = sb_religion, fill = sb_religion)) +
  geom_smooth(method = "lm", se = TRUE, alpha = 0.3) +
  ylab("Predicted Opposition to Persecution (1-7)") + 
  xlab("Attitudes Towards China") +   
  scale_y_continuous(limits = c(4.05, 5.3), breaks = seq(4.1, 5.3, by = 0.3)) +
  theme_bw() + 
  labs(title="", color="") + 
  theme(plot.title = element_text(size = 8, hjust = 0.5), 
        legend.position="bottom", 
        legend.title = element_blank()) +
  scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9", "#009E73")) +
  scale_fill_manual(values = c("gray","gray","gray","gray")) +
  guides(fill = "none")

china_moderation_plot

#ggsave("Figure4_Predicted_China_Therm.png",china_moderation_plot, width=6.5 ,height = 5, dpi=600)
#######################################################################################################

#######################################################################################################
#### Figure 5: Conditional Average Treatment Effects Across Christian Identity ####

data$christian <- ifelse(
  data$religious_affiliation %in% c(
    "Protestant or Evangelical", 
    "Roman catholic", 
    "Roman catholic,Protestant or Evangelical", 
    "Protestant or Evangelical,Orthodox", 
    "Orthodox"
  ), 
  1, 0
)

table(data$christian)

model_re_moderation <- lmer(index ~ as.factor(sb_religion)*as.factor(christian) + (1 | Country), data = data)
summary(model_re_moderation)

model_fe_moderation <- lm(index ~  as.factor(sb_religion)*as.factor(christian) + as.factor(Country), data=data) 
summary(model_fe_moderation)

## Table A12: Christian Identity Moderation Models##
moderation_models <- list(model_re_moderation, model_fe_moderation)

modelsummary(moderation_models, output = "latex", coef_map = , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('$^a$' = .1, '*'=0.05))


## Figure 5: Predicted Opposition Across Respondent Christianity ##
predictions_data <- plot_predictions(model_fe_moderation, by=c("christian","sb_religion"), draw = FALSE)

christian_predicted_plot <- ggplot(predictions_data, aes(x = factor(christian), y = estimate, color = sb_religion)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), size = 0.5, linewidth = 0.5,position=position_dodge(width=0.25)) +
  ylab("Predicted Opposition to Persecution (1-7)") +
  scale_x_discrete(breaks = c(0, 1), labels = c("Not Christian", "Christian")) +
  scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9", "#009E73")) +
  theme_bw() +
  theme(plot.title = element_text(size = 8, hjust = 0.5),
        legend.position="bottom") +
  labs(color="", x = "")

christian_predicted_plot

#ggsave("Figure5_Predicted_Christian.png",christian_predicted_plot, width=6.5 ,height = 5, dpi=600)
#######################################################################################################

#######################################################################################################
#### Table A13: Regression Models for Individual Outcomes Robustness  ####

model_re_worried <- lmer(persecution_stop ~ as.factor(sb_religion) + (1 | Country), data = data)
model_fe_worried <- lm(persecution_stop ~ as.factor(sb_religion) + as.factor(Country), data = data)

model_re_stop <- lmer(ch_govt_stop ~ as.factor(sb_religion) + (1 | Country), data = data)
model_fe_stop <- lm(ch_govt_stop  ~ as.factor(sb_religion) + as.factor(Country), data = data)
summary(model_re_stop)

model_re_country <- lmer(govt_pressure ~ as.factor(sb_religion) + (1 | Country), data = data)
model_fe_country <- lm(govt_pressure ~ as.factor(sb_religion) + as.factor(Country), data = data)
summary(model_re_country)

model_re_UN <- lmer(un_pressure ~ as.factor(sb_religion) + (1 | Country), data = data)
model_fe_UN <- lm(un_pressure ~ as.factor(sb_religion) + as.factor(Country), data = data)
summary(model_re_UN)

model_re_behavior <- lmer(petition ~ as.factor(sb_religion) + (1 | Country), data = data)
model_fe_behavior <- lm(petition ~ as.factor(sb_religion) + as.factor(Country), data = data)
summary(model_fe_behavior)


models <- list(model_re_worried,model_fe_worried, model_re_stop, model_fe_stop,model_re_country,model_fe_country,model_re_UN,model_fe_UN,model_re_behavior,model_fe_behavior)
modelsummary(models, output = "latex", coef_map = , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))
#######################################################################################################

#######################################################################################################
#### Table A14: Random Slopes and Random Intercept Models #### 

model_re_index <- lme(index ~ sb_religion, 
                      random = ~ sb_religion | Country, 
                      data = data,
                      method = "REML",
                      control = lmeControl(opt = "optim"))

model_re_worried <- lme(persecution_stop ~ sb_religion, 
                        random = ~ sb_religion | Country, 
                        data = data,
                        method = "REML",
                        control = lmeControl(opt = "optim"))

model_re_stop <- lme(ch_govt_stop ~ sb_religion, 
                     random = ~ sb_religion | Country, 
                     data = data,
                     method = "REML",
                     control = lmeControl(opt = "optim"))

model_re_country  <- lme(govt_pressure ~ sb_religion, 
                         random = ~ sb_religion | Country, 
                         data = data,
                         method = "REML",
                         control = lmeControl(opt = "optim"))

model_re_UN  <- lme(un_pressure ~ sb_religion, 
                    random = ~ sb_religion | Country, 
                    data = data,
                    method = "REML",
                    control = lmeControl(opt = "optim"))

model_re_behavior  <- lme(petition ~ sb_religion, 
                          random = ~ sb_religion | Country, 
                          data = data,
                          method = "REML",
                          control = lmeControl(opt = "optim"))


models <- list(model_re_index, model_re_worried, model_re_stop, model_re_country,model_re_UN,model_re_behavior)

modelsummary(models, output = "latex", coef_map = , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))
#######################################################################################################

#######################################################################################################
#### Table A15: Fixed Effects Models (HC2 Standard Errors) ####

model_fe_robust <- lm_robust(index ~ as.factor(sb_religion) + as.factor(Country), data = data)
summary(model_fe_robust)

model_fe_worried_robust <- lm_robust(persecution_stop ~ as.factor(sb_religion) + as.factor(Country), data = data)
summary(model_fe_worried_robust)

model_fe_stop_robust <- lm_robust(ch_govt_stop  ~ as.factor(sb_religion) + as.factor(Country), data = data)
summary(model_fe_stop_robust)

model_fe_country_robust <- lm_robust(govt_pressure ~ as.factor(sb_religion) + as.factor(Country), data = data)
summary(model_fe_country_robust)

model_fe_UN_robust <- lm_robust(un_pressure ~ as.factor(sb_religion) + as.factor(Country), data = data)
summary(model_fe_UN_robust)

model_fe_behavior_robust <- lm_robust(petition ~ as.factor(sb_religion) + as.factor(Country), data = data)
summary(model_fe_behavior_robust)

fe_robust_models <- list(model_fe_robust,model_fe_worried_robust, model_fe_stop_robust, model_fe_country_robust,model_fe_UN_robust,model_fe_behavior_robust)

modelsummary(fe_robust_models, output = "latex", coef_map = , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))
#######################################################################################################

#######################################################################################################
#### Table A16: Random Intercept Models (CR2 Standard Errors) ####

# WARNING: Code takes a while to run...# 

model_re_robust <- lmer(index ~ as.factor(sb_religion) + (1 | Country), data = data)
robust_se_1 <- coef_test(model_re_robust, vcov = "CR2", cluster=data$Country)
print(robust_se_1)

model_re_worried_robust <- lmer(persecution_stop ~ as.factor(sb_religion) + (1 | Country), data = data)
robust_se_2<- coef_test(model_re_worried_robust, vcov = "CR2", cluster=data$Country)
print(robust_se_2)

model_re_stop_robust <- lmer(ch_govt_stop ~ as.factor(sb_religion) + (1 | Country), data = data)
robust_se_3<- coef_test(model_re_stop_robust, vcov = "CR2", cluster=data$Country)
print(robust_se_3)

model_re_country_robust <- lmer(govt_pressure ~ as.factor(sb_religion) + (1 | Country), data = data)
robust_se_4 <- coef_test(model_re_country_robust, vcov = "CR2", cluster=data$Country)
print(robust_se_4)

model_re_UN_robust <- lmer(un_pressure ~ as.factor(sb_religion) + (1 | Country), data = data)
robust_se_5 <- coef_test(model_re_UN_robust, vcov = "CR2", cluster=data$Country)
print(robust_se_5)

model_re_behavior_robust <- lmer(petition ~ as.factor(sb_religion) + (1 | Country), data = data)
robust_se_6 <- coef_test(model_re_behavior_robust, vcov = "CR2", cluster=data$Country)
print(robust_se_6)

re_robust_models <- list(model_re_robust,model_re_worried_robust, model_re_stop_robust, model_re_country_robust, model_re_UN_robust,model_re_behavior_robust)
robust_se_list <- list(robust_se_1, robust_se_2, robust_se_3,  robust_se_4, robust_se_5, robust_se_6)

robust_se_1_SE <- robust_se_1[, "SE"]
robust_se_1_pvalues <- robust_se_1[, "p_Satt"]

robust_se_2_SE <- robust_se_2[, "SE"]
robust_se_2_pvalues <- robust_se_2[, "p_Satt"]

robust_se_3_SE <- robust_se_3[, "SE"]
robust_se_3_pvalues <- robust_se_3[, "p_Satt"]

robust_se_4_SE <- robust_se_4[, "SE"]
robust_se_4_pvalues <- robust_se_4[, "p_Satt"]

robust_se_5_SE <- robust_se_5[, "SE"]
robust_se_5_pvalues <- robust_se_5[, "p_Satt"]

robust_se_6_SE <- robust_se_6[, "SE"]
robust_se_6_pvalues <- robust_se_6[, "p_Satt"]

texreg(
  list(model_re_robust, model_re_worried_robust, model_re_stop_robust, 
       model_re_country_robust, model_re_UN_robust, model_re_behavior_robust),
  override.se = list(robust_se_1_SE, robust_se_2_SE, robust_se_3_SE, 
                     robust_se_4_SE, robust_se_5_SE, robust_se_6_SE),
  override.pvalues = list(robust_se_1_pvalues, robust_se_2_pvalues, robust_se_3_pvalues, 
                          robust_se_4_pvalues, robust_se_5_pvalues, robust_se_6_pvalues),
  stars = c(0.05), dcolumn = TRUE, booktabs = TRUE, digits = 3
)
#######################################################################################################

#######################################################################################################
#### Conditional Average Treatment Effects Across Other Foreign Policy Views on China #### 

model_re_moderation_economy <- lmer(index ~ as.factor(sb_religion)*China_important_economy + (1 | Country), data = data)
summary(model_re_moderation)

model_fe_moderation_economy <- lm(index ~ sb_religion*China_important_economy + as.factor(Country), data = data)
summary(model_fe_moderation)

prediction_data <- plot_predictions(model_fe_moderation_economy, by=c("China_important_economy","sb_religion"), draw = FALSE)

library(ggplot2)
china_important_economy_plot <- ggplot(prediction_data, aes(x = China_important_economy, y = estimate, color = sb_religion, fill = sb_religion)) +
  geom_smooth(method = "lm", se = TRUE, alpha = 0.3) +
  ylab("") + 
  xlab("China Important for Country Economy") +   
  scale_y_continuous(limits = c(3.1, 5.65), breaks = seq(3.1, 5.65, by = 0.4)) +
  scale_x_continuous(limits = c(1,7), breaks=seq(1,7,by=1)) +
  theme_bw() + 
  labs(title="", color="") + 
  theme(plot.title = element_text(size = 8, hjust = 0.5), 
        legend.position="bottom", 
        legend.title = element_blank()) +
  scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9", "#009E73")) +
  scale_fill_manual(values = c("gray","gray","gray","gray")) +
  guides(fill = "none")

model_re_moderation_foreign_policy <- lmer(index ~ as.factor(sb_religion)*China_align_foreign_policy + (1 | Country), data = data)
summary(model_re_moderation)

model_fe_moderation_foreign_policy <- lm(index ~ sb_religion*China_align_foreign_policy + as.factor(Country), data = data)
summary(model_fe_moderation)

prediction_data <- plot_predictions(model_fe_moderation_foreign_policy, by=c("China_align_foreign_policy","sb_religion"), draw = FALSE)

align_foreign_policy_plot <- ggplot(prediction_data, aes(x = China_align_foreign_policy, y = estimate, color = sb_religion, fill = sb_religion)) +
  geom_smooth(method = "lm", se = TRUE, alpha = 0.3) +
  ylab("Predicted Opposition to Persecution") + 
  xlab("Country Foreign Policy Align with China") +   
  scale_y_continuous(limits = c(3.1, 5.65), breaks = seq(3.1, 5.65, by = 0.4)) +
  scale_x_continuous(limits = c(0,10), breaks=seq(0,10,by=1)) +
  theme_bw() + 
  labs(title="", color="") + 
  theme(plot.title = element_text(size = 8, hjust = 0.5), 
        legend.position="bottom", 
        legend.title = element_blank()) +
  scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9", "#009E73")) +
  scale_fill_manual(values = c("gray","gray","gray","gray")) +
  guides(fill = "none")

align_foreign_policy_plot

model_re_moderation_human_rights <- lmer(index ~ as.factor(sb_religion)*country_advance_human_rights_China + (1 | Country), data = data)
summary(model_re_moderation)

model_fe_moderation_human_rights <- lm(index ~ sb_religion*country_advance_human_rights_China + as.factor(Country), data = data)
summary(model_fe_moderation)

prediction_data <- plot_predictions(model_fe_moderation_human_rights, by=c("country_advance_human_rights_China","sb_religion"), draw = FALSE)

country_advance_hr_plot <- ggplot(prediction_data, aes(x = country_advance_human_rights_China, y = estimate, color = sb_religion, fill = sb_religion)) +
  geom_smooth(method = "lm", se = TRUE, alpha = 0.3) +
  ylab("Predicted Opposition to Persecution") + 
  xlab("Country Advance Human Rights in China") +   
  scale_y_continuous(limits = c(3.1, 5.65), breaks = seq(3.1, 5.65, by = 0.4)) +
  scale_x_continuous(limits = c(1,7), breaks = seq(1,7)) +
  theme_bw() + 
  labs(title="", color="") + 
  theme(plot.title = element_text(size = 8, hjust = 0.5), 
        legend.position="bottom", 
        legend.title = element_blank()) +
  scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9", "#009E73")) +
  scale_fill_manual(values = c("gray","gray","gray","gray")) +
  guides(fill = "none")

model_re_moderation_major_power <- lmer(index ~ as.factor(sb_religion)*feeling_China_major_power + (1 | Country), data = data)
summary(model_re_moderation)

model_fe_moderation_major_power  <- lm(index ~ sb_religion*feeling_China_major_power + as.factor(Country), data = data)
summary(model_fe_moderation)

prediction_data <- plot_predictions(model_fe_moderation_major_power , by=c("feeling_China_major_power","sb_religion"), draw = FALSE)

china_great_power_plot <- ggplot(prediction_data, aes(x = feeling_China_major_power, y = estimate, color = sb_religion, fill = sb_religion)) +
  geom_smooth(method = "lm", se = TRUE, alpha = 0.3) +
  ylab("") + 
  xlab("Attitudes Towards China's Rise as Major Power") +   
  scale_y_continuous(limits = c(3.1, 5.65), breaks = seq(3.1, 5.65, by = 0.4)) +
  theme_bw() + 
  labs(title="", color="") + 
  theme(plot.title = element_text(size = 8, hjust = 0.5), 
        legend.position="bottom", 
        legend.title = element_blank()) +
  scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9", "#009E73")) +
  scale_fill_manual(values = c("gray","gray","gray","gray")) +
  guides(fill = "none")

combined <- ggarrange(align_foreign_policy_plot, china_important_economy_plot, country_advance_hr_plot, china_great_power_plot, ncol = 2, nrow = 2, common.legend = TRUE, legend = "bottom")  

#ggsave("FigureA1_Additional_China_Moderations.jpeg",combined, width=8 ,height = 6, dpi=1000)

## Tables A17 and A18-19 ##
re_models <- list(model_re_moderation_economy,model_re_moderation_foreign_policy, model_re_moderation_human_rights, model_re_moderation_major_power)

fe_models <- list(model_fe_moderation_economy,model_fe_moderation_foreign_policy, model_fe_moderation_human_rights, model_fe_moderation_major_power)

modelsummary(re_models, output = "latex", coef_map = , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))

modelsummary(fe_models, output = "latex", coef_map = , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))
#######################################################################################################

#######################################################################################################
#### Conditional Average Treatment Effects Across National and European Identity ####

nat_identity <- select(data, nationality_important_identity, glad_to_be_nationality, bond_people_country)

alpha(nat_identity)

data$national_identity <- (data$nationality_important_identity+data$glad_to_be_nationality+data$bond_people_country)/3 

european_identity <- select(data, European_important_identity, glad_to_be_European, bond_people_European)

alpha(european_identity)

data$european_identity <- (data$European_important_identity+data$glad_to_be_European+data$bond_people_European)/3

model_re_moderation_ni <- lmer(index ~ as.factor(sb_religion)*national_identity + (1 | Country), data = data)
summary(model_re_moderation_ni)

model_fe_moderation_ni <- lm(index ~ as.factor(sb_religion)*national_identity + as.factor(Country), data = data)
summary(model_fe_moderation_ni)

model_re_moderation_ei <- lmer(index ~ as.factor(sb_religion)*european_identity + (1 | Country), data = data)
summary(model_re_moderation_ei)

model_fe_moderation_ei <- lm(index ~ as.factor(sb_religion)*european_identity +as.factor(Country), data = data)
summary(model_fe_moderation_ei)

mean_nat_ident <- mean(data$national_identity)
sd_nat_ident <- sd(data$national_identity)

mean_euro_ident <- mean(data$european_identity)
sd_euro_ident <- sd(data$european_identity)

## Table A20: National and European Identity Moderation Models##
identity_moderation_models <- list(model_re_moderation_ni,model_fe_moderation_ni,model_re_moderation_ei,model_fe_moderation_ei)

modelsummary(identity_moderation_models, output = "latex", coef_map = , gof_omit = "AIC|BIC|Log.Lik.|RMSE|F",stars = c('*' = .05))


## Figure A2: Predicted Opposition Across National and European Identity Attachments ##
prediction_data_ni <- avg_predictions(model_fe_moderation_ni, variables=list(national_identity=c(mean_nat_ident-sd_nat_ident,mean_nat_ident,mean_nat_ident+sd_nat_ident)), by=c("sb_religion","national_identity"))


prediction_data_ei <- avg_predictions(model_fe_moderation_ei, variables=list(european_identity=c(mean_euro_ident-sd_euro_ident,mean_euro_ident,mean_euro_ident+sd_euro_ident)), by=c("sb_religion","european_identity"))

national_identity_predicted <- ggplot(prediction_data_ni, 
                                      aes(x = factor(national_identity), 
                                          y = estimate, 
                                          color = sb_religion)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), size = 0.4, linewidth = 0.5, position=position_dodge(width=0.25)) +
  ylab("Predicted Opposition to Persecution") +
  xlab("National Identity Attachment") +
  scale_x_discrete(labels = c("-1sd", "Mean", "+1sd")) +
  scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9", "#009E73")) +
  scale_y_continuous(limits = c(4.3, 5.2), breaks = seq(4.3, 5.2, by = 0.2)) +
  theme_bw() +
  theme(plot.title = element_text(size = 8, hjust = 0.5),
        legend.position = "bottom",
        legend.title = element_blank()) +
  labs(color = "")

national_identity_predicted

cor(data$national_identity,data$european_identity)

european_identity_predicted <- ggplot(prediction_data_ei, 
                                      aes(x = factor(european_identity), 
                                          y = estimate, 
                                          color = sb_religion)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), 
                  size = 0.4, linewidth = 0.5, position = position_dodge(width = 0.25)) +
  ylab("") +  
  xlab("European Identity Attachment") +
  scale_x_discrete(labels = c("-1sd", "Mean", "+1sd")) +
  scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9", "#009E73")) +
  scale_y_continuous(limits = c(4.3, 5.2), breaks = seq(4.3, 5.2, by = 0.2)) +
  theme_bw() +
  theme(plot.title = element_text(size = 8, hjust = 0.5),
        legend.position = "bottom",
        legend.title = element_blank()) +
  labs(color = "")

european_identity_predicted

identity_plot <- ggarrange(national_identity_predicted, european_identity_predicted, nrow=1, common.legend = TRUE, legend="bottom")

identity_plot

#ggsave("FigureA2_Predicted_Identity_Attachment.jpeg",identity_plot, width=6.5 ,height = 5, dpi=1000)
#######################################################################################################

#######################################################################################################
#### Descriptives and Balance Tables  ####

split_datasets <- split(data, data$Country)
list2env(split_datasets, envir = .GlobalEnv)

balance_summaries <- lapply(split_datasets, function(data) {
  datasummary_balance(Gender + Age + education + religious_affiliation + religious_attendance ~ sb_religion, 
                      dinm = TRUE, dinm_statistic = "p.value", 
                      data = data, output = "latex")
})

names(balance_summaries) <- names(split_datasets)

balance_summaries[["CZ"]]
balance_summaries[["FR"]]
balance_summaries[["GE"]]
balance_summaries[["HU"]]
balance_summaries[["IT"]]
balance_summaries[["LV"]]
balance_summaries[["PL"]]
balance_summaries[["RU"]]
balance_summaries[["SR"]]
balance_summaries[["SK"]]
balance_summaries[["ES"]]
balance_summaries[["SW"]]
balance_summaries[["UK"]]
#######################################################################################################